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ABSTRACT 

The disruption of stars by supermassive black holes has been linked to more than a dozen flares in the cores 
of galaxies out to redshift z ~ 0.4. Modeling these flares properly requires a prediction of the rate of mass re- 
turn to the black hole after a disruption. Through hydrodynamical simulation, we show that aside from the full 
disruption of a solar mass star at the exact limit where the star is destroyed, the common assumptions used to 
estimate M(f), the rate of mass return to the black hole, are largely invalid. While the analytical approximation 
to tidal disruption predicts that the least-centrally concentrated stars and the deepest encounters should have 
more quickly-peaked flares, we find that the most-centrally concentrated stars have the quickest-peaking flares, 
and the trend between the time of peak and the impact parameter for deeply-penetrating encounters reverses 
beyond the critical distance at which the star is completely destroyed. We also show that the most-centrally 
concentrated stars produced a characteristic drop in M(t) shortly after peak when a star is only partially dis- 
rupted, with the power law index n being as extreme as -4 in the months immediately following the peak of a 
flare. Additionally, we find that n asymptotes to ~ -2.2 for both low- and high-mass stars for approximately 
half of all stellar disruptions. Both of these results are significantly steeper than the typically assumed n = -5/3. 
As these precipitous decay rates are only seen for events in which a stellar core survives the disruption, they 
can be used to determine if an observed tidal disruption flare produced a surviving remnant. We provide fitting 
formulae for four fundamental quantities of tidal disruption as functions of the star's distance to the black hole 
at pericenter and its stellar structure: The total mass lost, the time of peak, the accretion rate at peak, and the 
power-law index shortly after peak. These results should be taken into consideration when flares arising from 
tidal disruptions are modeled. 

Subject headings: accretion, accretion disks — black hole physics — gravitation — hydrodynamics — meth- 
ods: numerical 



1. INTRODUCTION 

Supermassive black holes (SMBHs) have been found to re- 
side at the centers of most galaxies. These black holes are or- 
bited by a cluster of stars that interact with one another gravi- 
tationally through stochastic encounters. Occasionally, an en- 
counter will shift a star onto an orbit that takes it within its 
tidal radius, defined as the distance at which the black hole's 
tidal forces would overcome the star's self-gravity at its sur- 
face (Frank & Rees 1976). A fraction of the star's mass then 
becomes bound to the black hole, and proceeds to fall back 
towards the star's original pericenter, eventually forming an 
accretion disk that results in a luminous flare with a luminos- 
ity comparable to the Eddington luminosity. 

The standard model of tidal disruption presumes that the 
star is completely destroyed, resulting in approximately half 
of the star's original mass falling back onto the black hole, 
with the debris possessing a variety of orbital periods result- 
ing from a spread of orbital energy that is "frozen in" at peri- 
center. First described in Rees (1988), the rate of fallback 
has been estimated both through increasingly sophisticated 
numerical simulations and analytical models. While previ- 
ous results have provided reasonable models for the fallback 
resulting from the complete disruptions of stars at the tidal ra- 
dius r t = /^(Mh/M*) 1 / 3 , where and are the mass and 
radius of the star and Mh is the mass of the black hole, they 
completely neglect partial stellar disruptions, in which a stel- 
lar core survives the encounter and only a fraction of the star's 
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mass becomes immediately bound to the black hole. These 
events are likely to be much more common than their com- 
plete disruption counterparts, both for the reason that the rate 
of encounters interior to the pericenter distance r p scales as 
r p (Hills 1988), and also that the disrupted star may return on 
subsequent orbits and be subject to disruption and/or further 
tidal dissipation. Additionally, many previous studies have fo- 
cused on stars of a single structural profile, usually selected to 
match the familiar profile of our own Sun. However, standard 
stellar mass functions predict that low-mass main sequence 
(MS) stars are more common (e.g. Kroupa et al. 1993), and 
thus may contribute significantly to the overall disruption rate. 
These stars are significantly less centrally concentrated than 
their solar mass brethren. 

The dynamics of stellar tidal disruption have been mod- 
eled by many authors using simple analytical arguments (Rees 
1988; Phinney 1989; Lodato et al. 2009; Kasen & Ramirez- 
Ruiz 2010), increasingly complex dynamical models (Lu- 
minet & Marck 1985; Carter & Luminet 1985; Luminet & 
Carter 1986; Diener et al. 1995; Ivanov & Novikov 2001), and 
hydrodymical simulations utilizing either an Eulerian (Evans 
& Kochanek 1989; Khokhlov et al. 1993a,b; Diener et al. 
1997; Guillochon et al. 2009) or Lagrangian (Nolthenius & 
Katz 1982; Bicknell & Gingold 1983; Laguna et al. 1993; 
Kobayashi et al. 2004; Rosswog et al. 2008; Lodato et al. 
2009; Ramirez-Ruiz & Rosswog 2009; Rosswog et al. 2009; 
Antonini et al. 2011) approaches. Very few of these studies 
have presented the effect varying r p on the amount of mass 

lost by the star, AM, or the effect on M(f), the rate at which 
the mass liberated from the star returns to pericenter. Given 
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that the viscous time is expected to be significantly shorter 
than the period of the returning debris, this M(t) is expected 
to track the luminosity L(t) closely. As the number of ob- 
served disruptions increases, and as the cadence and quality 
of data improves, it becomes increasingly more important to 
improve models of M(t) for disruptions of all kinds. 

In this paper, we present the results of 43 hydrodynamical 
simulations at high-resolution representing the disruption of 
both low-mass and high-mass MS stars. This provides, for 
the first time, a complete picture of the feeding of SMBHs 
by the disruption of MS stars. While the expected trend of 
smaller mass accretion rates for progressively more grazing 
encounters is reproduced, our study reveals several surprises 
on how disruptions work, particularly on the effect of stellar 
structure and how the fallback rate scales for both grazing and 
deep encounters. Contrary to what is expected from the freez- 
ing model, in which only the distribution of mass at pericenter 
is considered, the non-linear response of the star to the tidal 
field is found to play a crucial role in determining M(t). Our 
simulations show that the simple models previously employed 
to predict the rate of fallback do not capture the full dynamics 
of the problem, and are only appropriate for anything other 
than the full disruption at exactly the tidal radius. 

We find that the decay rate of M{t) does not settle to a con- 
stant value until a few months after the disruption for all dis- 
ruptions by black holes with Mh > 10 6 M Q , implying that the 
range of characteristic decay rates used to identify tidal dis- 
ruption flares should be widened to include events that may 
not follow the fiducial f" 5 / 3 decay rate. For partial disrup- 
tions, the decay rate at a few years after the disruption de- 
pends crucially on the hydrodynamical evolution of the de- 
bris stream. This means that simulations must cover more 
than a few stellar dynamical timescales after the disruption, 
with the final functional form of M{t) not being established 
until the star is many hundreds of tidal radii away from the 
black hole. And while we do find that there are differences be- 
tween the fallback functions calculated for the disruptions of 
profiles characteristic of low- and high-mass stars, the mass- 
radius relationship of MS stars results in a family of fall- 
back curves that are difficult to distinguish from one another 
for stars of O.3M > M* > 1.0M Q without considering sec- 
ondary features related to the shape of the fallback curves 
themselves, such as the decay rate of M(t), characterized by a 
time-dependent power-law index n(f). 

This paper is organized as follows. We describe our numer- 
ical method, initial models, and method for calculating AM 
and M{t) in Section 2. The results of these simulations and 
how they improve our understanding of stellar tidal disrup- 
tions is described in Section 3. A discussion of the general 
trends and their effect on the observable features of tidal dis- 
ruptions is presented in Section 4. Finally, we provide fitting 
formula to four characteristic variables describing disruptions 
of stellar profiles characteristic of low- and high-mass stars in 
Appendix A. 

2. METHOD 

Our simulations of tidal disruption are performed in FLASH 
(Fryxell et al. 2000), an adaptive-mesh grid-based hydrody- 
namics code which includes self-gravity. The standard hydro- 
dynamical equations are solved using the directionally split 
piecewise-parabolic approach (Colella & Woodward 1984) as 
provided by the FLASH software, which has a small numeri- 
cal viscosity and diffusivity as compared to the available un- 



split solvers Therefore, the principle source of entropy gen- 
eration is via shocks, if they are present. The solution to the 
Riemann problem is sensitive to the velocity of the frame in 
which the problem is solved, and poor solutions can be re- 
turned in regions where the velocity of the frame is many 
times larger than the sound speed (Tasker et al. 2008; Robert- 
son et al. 2010; Springel 2010). As stars that are disrupted by 
SMBHs are traveling at a fraction of the speed of light c, we 
perform our simulations in the rest-frame of the star where the 
velocities are ~ ^/c s v p , where c s is the sound speed within the 
star and v p is the velocity at pericenter. 

Our method is very similar to what is presented in Guil- 
lochon et al. (2011), except that we now utilize version 4.0 
of the FLASH software, which has a greatly improved multi- 
pole gravity solver 2 . A key parameter of the multipole gravity 
solver is the maximum angular number of the multipole ex- 
pansion l m . A test simulation setting l m = 10 showed multi- 
ple recollapse points for a nearly-complete disruption, which 
is not expected to occur in disrupted stars (as described in 
Section 3.1). We suspected this behavior was a consequence 
of the large aspect ratio of the debris stream, which results 
in gravity being under-resolved if the multipole expansion is 
truncated at small l m . With l m set to 20, only a single rec- 
ollapse occurs, as is expected. For l m = 40, we found no 
significantly difference in any quantities of interest as com- 
pared to l m = 20, except for cases in which a very small rem- 
nant survives the disruption. In these cases, the error is in 
the mass of the surviving object, which is difficult to resolve 
for marginally surviving stars (see Section 3.1). The results 
presented in this work all use l m = 20 for optimal runtime ef- 
ficiency. 

Additionally, we add a truncation density parameter pdamp 
that is set to 10~ 18 g cm" 1 , a factor of 10 larger than the fluff 
density. Material with density less than this value is not in- 
cluded in the multipole expansion, nor are any gravitational 
forces applied to this material. This is necessary as the do- 
main is extremely large as compared to the initial star, and so 
even 1O _11 M of material can introduce a significant error in 
the calculated position of the center of mass. We also only 
include material with a density greater than 10% of the star's 
original peak density when calculating the location to use as 
the multipole expansion point; however all material with den- 
sity greater than pdamp is included when calculating the magni- 
tude of each of the multipole terms. This is done because the 
total center of mass does not always spatially coincide with 
the peak density, which can result in a multipole expansion 
that is a poor representation of the underlying density distri- 
bution. 

2.1. Parameter Study 

Ignoring general relativistic effects and stellar rotation, it 
may seem that a complete study of tidal disruptions would re- 
quire an exhaustive study of the various combinations of six 
parameters: M*, Mh, the orbital eccentricity <?, the poly- 
tropic index 7, and (3 = r t /r p . As an exhaustive search of a six- 
dimensional parameter space is prohibitive, we wish to reduce 
the number of free parameters to a more manageable number. 

1 II 

For fixed /3, both r p and v p scale as M h , and thus the peri- 
center crossing time f p is independent of M^. Additionally, as 
the mass ratio approaches infinity, the asymmetry of the tidal 
field becomes progressively less important as <C r t , with 

2 See http://flash.uchicago.edu for details 



FIG. 1 . — Snapshots of the density log p for all 7 = 4/3 simulations at t = 4 X 10 4 s after the start of each simulation, with white corresponding to the maximum 
density and black corresponding to 1CT 6 g cirT 1 . Each snapshot is labeled with the ratio of the tidal radius to the pericenter distance 0, The white arrows indicate 
the angle of the velocity vector at the time of each snapshot. The white dashed line separates simulations in which a core survives the encounter; although not 
visible here, recollapse does occur for the /3 = 1 .7 and 1 .8 simulations, but for I > 4x 10 4 s. 



the difference in the strength of the tidal field at pericenter 
between the near-side and the far-side for a 10 6 : 1 encounter 
being ~ 3% (Guillochon et al. 2011). And as most of the 
stars that are scattered into disruptive orbits originate from the 
sphere of influence or beyond (Magorrian & Tremaine 1999; 
Wang & Merritt 2004), the orbital eccentricity of almost all 
disrupted stars is approximately unity 3 . 

It then follows that as the ratio of th e time of t he encounter 
f p to the star's dynamical time fd yn = \/Rl/ GM* , shape of the 
orbit (which depends on e), and asymmetry of the tides are 
nearly identical for all encounters of interest for a fixed /?, the 
tidal force applied to the star as a function of time is indepen- 
dent of Mh, <?, , and . Thus, we find that the vast majority 
of stellar disruptions by SMBHs can be described by just two 
parameters: f3 and 7, with all other parameters obeying sim- 
ple scaling relations. While previous numerical studies have 
considered the effect of varying 7 on M(t )(Rosswog et al. 
2009; Lodato et al. 2009; Ramirez-Ruiz & Rosswog 2009), 
the present work is the first to explore the effect of varying /? 
on M(t) in cases ranging from no mass loss to deeply pene- 
trating encounters 4 . 

To explore this reduced but physically motivated parameter 
space, we run a series of simulations assuming = M Q and 
Mh = 10 6 M Q . Our stars are constructed as polytropes, with 
the polytropic 7 being set to either 5/3 or 4/3, representative 
of both low- and high-mass stars, respectively. During the 
simulation, the stars are evolved hydrodynamically accord- 
ing to a T = 5/3 equation of state, with the difference be- 

3 See (Madigan et al. 201 1) for a discussion of resonant relation processes 
that may produce a different distribution of eccentricities for stellar disrup- 
tions. 

4 Note that Laguna et al. (1993) do present M(t) from low-resolution sim- 
ulations for three different /3 values, two of which are very deeply penetrating 

(fi > 5). 



tween 7 and T for high-mass stars being a consequence of 
the radiation transfer within the star (Chandrasekhar 1939). 
These one-dimensional profiles are then mapped to the three- 
dimensional grid, with initially uniform refinement across the 
star. The star is then relaxed for 10 4 s at the center of a cu- 
bical domain, which is 4 x 10 14 cm on a side. The domain 
is initially composed of a single 8 3 block, which is then bi- 
sected into smaller 8 3 blocks as many as 15 times, resulting in 
a minimum cell size of 3 x 10 9 cm, or approximately 2% of 
the star's original diameter. Our refinement criteria is solely 
dependent on the density relative to the initial central density, 
with a factor of two reduction in resolution for each factor of 
a hundred in density. Regions within the simulation that are 
within 1 % of the peak density are always maximally refined. 

We ran simulations for 23 different impact parameters (3 = 
r t /r p ranging from 0.6 to 4.0 for 7 = 4/3, and 20 different f3 
ranging from 0.5 to 2.5 for 7 = 5/3. Two additional simula- 
tions were run at /3 = 0.5 for 7 = 4/3 and f3 = 0.45 for 7 = 5/3; 
as less than 10~ 6 M Q is observed to be removed from the stars 
in these two borderline cases, we conlude that no mass is lost 
for values of (3 less than the above quoted ranges in /3. Snap- 
shots from each simulation recorded shortly after pericenter 
are shown in Figures 1 and 2. 

2.2. Calculation of AM and M{t) 

Our hydrodynamical simulations enable us to calculate the 
binding energy of the material to the black hole dM/dE. This 
function can be used to determine the feeding rate as a func- 
tion of time through Kepler's third law, 

dMdE 
dE dt 



M(t): 



27r „ dM <;/, 

— (GM h ) 2 / 3 — r 5 / 3 . 

3 dE 



(1) 



For full disruptions, the entirety of the star's original mass 
is included in the calculation of dM/dE, approximately half 
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FIG. 2. — Snapshots of log p for all 7 = 5/3 simulations at t = 4 X 10 4 s, colors 
does occur for the = 0.75 - 0.85 simulations, but for (>4x 10 4 s. 

of which will have specific orbit energy E > and is thus un- 
bound from the black hole. For partial disruptions, the criteria 
for determining which material to include in the determination 
of dM/dE is less straightforward, as what will be accreted by 
the black hole is only the material that the star's gravity is un- 
able to retain. As the star is on a parabolic orbit, the distance 
from the black hole changes rapidly as a function of time, and 

thus the star's Hill radius «h(0 = KMb OU nd(0 / Mb) 1 ' 3 is also 
time-dependent, introducing some ambiguity into the deter- 
mination of the self-bound mass Mb 0U nd(?)- 

In principle, the distance of matter from the surviving star 
can be compared to «h(0 to determine what mass is bound to 
the star. However, the continual reaccretion of matter means 
that the star is extended, non-spherical, and dynamically unre- 
laxed for many dynamical timescales, and thus the appropriate 
mass to use in the calculation of oh is uncertain. To circum- 
vent this, we choose an iterative energy-based approach that 
we find converges quickly to a solution. First, we calculate 
the material that remains bound to the star, where the initial 
reference point is taken to be at the location of the star's peak 
density, which has a velocity v pea k- The specific binding en- 
ergy of material in a given cell is calculated as 

1 2 

E*,i= - (v, -v peak ) (2) 

where is the gravitational self -potential as calculated by the 
multipole solver. The center of momentum v cm is then deter- 
mined by summing over all mass elements for which E* i < 

22s^<aPi v t 

where p, and V; are the cell density and volume. Equation (2) 
is then re-evaluated with v pea k being replaced by v cm . This 
process is repeated until v cm (and thus M(, oun d) converges to a 
constant value. While this approach yields a value for Mb OU nd 
in most cases, the question of whether an object is completely 
destroyed is somewhat complicated by the fact that the tidal 
force formally approaches zero as Ar — > 0, and thus there is 
always some material for which v,- = v cm , resulting in a in- 
finitesimal, but non-zero Mbound even as the time since disrup- 
tion t—t& — y 00, 

Figure 3 shows how the maximum density within a sim- 
ulation p max compares to the star's initial maximum density 
Pmax.o for six simulations (three for 7 = 4/3 and three for 



and annotations are the same as in Figure 1 . Although not visible here, recollapse 

7 = 5/3). Two of the simulations shown for each 7 exhibit 
a continuous decrease in p max , showing no signs of recol- 
lapse, whereas the third simulation for each 7 shows an in- 
crease in density sometime after pericenter, eventually set- 
tling to a constant value as the collapsed object dynamically 
relaxes. As a check on the convergence of Mb OU nd for all the 
simulations presented in this work, we compute the quantity 

T = I Mbound /Abound I (t — fd), which expresses the fractional 
change in Mb oun d since the time of disruption. Disruptions in 
which a self-bound core forms asymptote quickly to a con- 
stant Mbound, and thus small values of T, whereas disruptions 
in which p max consistently decreases show T ~ 1 at all times. 
The only disruptions in which the final core mass has not 
completely converged are the borderline survival cases (e.g. 
(3 slightly less than /3<j). However, while the fractional error 
in Mb olm d is large for the borderline cases, the definition of 
AM = M.J, - Mbound means that the amount of mass lost from 
the star (and also the amount of mass bound to the black hole) 
is well-determined. 

Once v cm has been determined, all material for which E*^ < 
is excluded, and the binding energy to the black hole E is 
calculated. This data is then binned in E, the result of which 
is used to determine dM/dE. The values of AM and M(t) 
presented in the figures in the subsequent sections are all gen- 
erated from snapshots that are produced at t = 2.5 x 10 5 s af- 
ter the start of each simulation (unless otherwise noted), or 
approximately 100 dynamical times after pericenter. 

3. HYDRODYNAMICS OF THE TIDAL DISRUPTION OF MS STARS 

Many assumptions about the way partial and full disrup- 
tions work have never been tested beyond analytical approxi- 
mations. Quantities that have been estimated include the time 
of return of the most bound material f most , the time of peak 
accretion rate t p&: ± and the magnitude of this rate Mp ea k, and 
the amount of mass bound both to the star and to the black 
hole after the encounter. Additionally, it has always been 
presumed that the late-time evolution of the fallback con- 
verges to the f~ 5 / 3 decay law, whereas this is not necessarily 
true in partial disruptions where the surviving core may affect 
the binding energy of this material. We empirically measure 
these quantities from our calculations of dM/dE, and find that 
while some of the commonly-held assumptions are reason- 
ably accurate, many are not. Most of these assumptions arise 
from how the problem was originally formulated, in which 
the star's self-gravity is viewed as inconsequential, and only 
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FIG. 3. — Evolution of maximum density p raax and bound mass Mb olln( j as a function of time since disruption. In the left panel, the evolution of the ratio of 
Pmax to the original maximum density p max .o is shown for six simulations (filled regions), three for 7 = 4/3 (orange, solid lines) with /3 = 0.85, 0.9, and 0.95, and 
three for 7 = 5/3 (light blue, dashed lines) with /3 = 1.8, 1.85, and 1.9. In the right panel, the parameter T = lAfbound/^boundlC — 'd) is shown for all simulations, 
demonstrating the convergence of the calculate M(, olln[ j for all the simulations presented in this work (see Section 2.2 for details). When the value of this quantity 
is close to unity, Mb olmc i is still changing by order unity over that timescale, indicating that the final bound mass cannot yet be determined at that t . The thick lines 
show simulations for which the star is considering to be destroyed after the encounter, whereas the thin lines show simulations for which a surviving core forms. 



the spread in binding energy across the star at pericenter is 
relevant in determining the features of the resulting M(t). We 
find that the star's self-gravity is critical in determining the 
resulting M(f), even for encounters with pericenters that are 
many times deeper than the tidal radius. 

3.1. The Boundary Between Survival or Destruction 

A collection of non-interacting particles in the presence of a 
point mass potential will all follow Keplerian orbits, provided 
that no outside force acts upon them. This means that once 
both the star's gravity and pressure become unimportant at a 
time close to the star's closest approach to a black hole, the 
position and velocity of each mass element can be recorded, 
and the future orbits of each part of the debris stream can be 
determined. It has been presumed that this condition is satis- 
fied at r t , the distance at which the tidal force is greater than 
the self-gravitational force at the object's surface. This as- 
sumption is flawed in that the tidal radius as classically de- 
fined does not denote the distance at which the tidal force 
dominates self-gravity for any point within the star, but rather 
only at its surface. The conditions necessary for a polytrope 
to lose mass due to the presence of an external tidal force have 
been previously determined in the context of the Roche prob- 
lem, which considers when the tidal force at the surface of an 
object exceeds the self-gravitational force in a circular orbit 
(Aizenman 1968; Chandrasekhar 1969). However, again, this 
limit only informs us as to when we expect the object to be- 
gin losing mass, and not the distance at which the object is 
completely destroyed. Additionally, the Roche limit is eval- 
uated under the assumption of hydrostatic equilibrium, and 
presumes that the orbital velocity is equal to that of a circular 
orbit v c , resulting in a different dynamical response than for 
parabolic encounters in which the pericenter velocity is \/2v c . 

The question of whether a star survives depends not on the 
ability of tidal forces to remove some mass, but on whether 
these forces are overwhelming enough to disrupt the star's 
densest regions. Furthermore, even if a star experiences a 
seemingly complete disruption, the star may be capable of 
recollapse into a self-bound object after the encounter under 
the proper conditions. It has been shown that gamma-law 
equations of state stiffer than V = 2 can result in the recol- 



lapse of material within expanding thin streams for infinites- 
imally small masses (Chandrasekhar 1961; Lee & Ramirez- 
Ruiz 2007). As stars are well-approximated by Y < 5/3 equa- 
tions of state, these instabilities are not expected to appear in 
stellar disruptions, and thus recollapse is not guaranteed for 
all 8. 

The affine model, as introduced in Carter & Luminet 
(1985), improved upon the initial estimates provided by the 
Roche approach by including the effects of the dynamical 
tide, but while this approach is able to evaluate the distance at 
which distortions become non-linear, it is not capable of deter- 
mining the actual distance at which disruptions occur. Later, 
Diener et al. (1995) extended the affine approach to calculate 
the critical impact parameter for full disruption finding 
/3d = L12 for 7 = 4/3 and (3 d = 0.67 for 7 = 5/3 polytropes, 
where /3<j is the critical impact parameter at which complete 
disruption ensues. More recently, the affine formalism was 
improved upon further by modeling the star as a nested set 
of ellipsoids, each of which respond dynamically to the exter- 
nal tidal field (Ivanov & Novikov 2001; Ivanov et al. 2003). 
While this model is the first analytical approach to provide 
estimates for AM, the simplifying assumptions made regard- 
ing the treatment of self-gravity, pressure, and geometry does 
not guarantee that the true solution can be recovered via this 
approach. 

In Figure 4 we show the amount of mass lost AM = M„ - 
-Abound (measured at the end of each simulation) as a function 
of /3forboth7 = 4/3 and7 = 5/3, with comparisons to Ivanov 
& Novikov (2001) shown in black. Remarkably, the model of 
Ivanov & Novikov comes quite close to predicting the critical 
B value as measured by these simulations, despite the assump- 
tions made, and is able to recover reasonable values for AM, 
although the scaling between AM and B is somewhat steeper 
than what is observed in the simulations. In particular, we 
find that stars can survive encounters for larger values of j3 
than the nested affine model predicts. We speculate that the 
method presented in Ivanov & Novikov could be extended to 
calculate M(f) if the time at which each ellipsoid becomes 
unbound were recorded, which given the low computational 
burden of this approach could be used to perform more exten- 
sive parameter space studies. 
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FIG. 4. — Fits to AM, with the fits to the 7 = 4/3 models being shown by 
the solid colored circles, and fits to the 7 = 5/3 models being shown by the 
open colored circles. Predictions of AM from Ivanov & Novikov (2001) for 
both 7 = 4/3 and 7 = 5/3 are represented by the black symbols/curves. The 
color coding matches that of Figures 1 and 2, with the impact parameters 
/3j beyond which stars are considered to be destroyed being denoted by the 
colored dot-dashed lines. 

We observe that while some stars appear initially to be com- 
pletely destroyed, with their cores being disrupted along with 
their envelopes (Figures 1 and 2), the debris stream can often 
recollapse many dynamical timescales after the encounter, re- 
sulting in a small yet self-bound remnant. The mass of the 
remnant that results is roughly equal to the amount of mass 
contained within a sphere centered at the recollapse point 
and with a radius equal to the cylindrical radius of the debris 
stream S. 

For collapsing gaseous cylinders, spurious condensations 
as the result of the accumulation of numerical error may de- 
velop if the Jeans length is not properly resolved (Truelove 
et al. 1997), with the source of that error being exacerbated by 
an inexact determination of the gravitational potential (Jiang 
et al. 2013). Truelove et al. found that no spurious gravi- 
tational collapse occurs if the ratio J of the grid scale to the 
Jeans length Aj = \J ttc 2 / Gp, where c s is the sound speed and 
p is the density, is always less than 0.25 in all grid cells at all 
times. In all of our simulations, the width of the debris stream 
is comparable to the star's initial size, and the resolution in the 
densest portion of the stream as it condenses is equal to the 
resolution used to resolve the original star (~ 50 grid cells). 
Therefore, in the case in which a recollapse marginally occurs 
(i.e. J/S ~ 1), J ~ 0.02, satisfying the Truelove criteria. 

For 7 = 4/3, we find that stars are destroyed for /? > /?d = 
1.85, i.e. no self-bound stellar remnant is produced. To ver- 
ify that we are adequately resolving the boundary between 
survival and destruction, we ran a single 7 = 4/3, /3 = 1 . 8 sim- 
ulation at double the linear resolution, and found a recollapse 
that results in a bound remnant of only a few percent of a so- 
lar mass, slightly smaller than what is found using our fiducial 
resolution. As the mass of the surviving star nears zero, the 
resolution requirements become progressively more restric- 
tive, as even slight changes in the cylindrical density profile 
or gravitational potential can alter the time of recollapse, and 
thus the final bound mass. For 7 = 5/3, we find that stars are 
destroyed for (3 > /3 d = 0.9. 

Numerical challenges aside, the exact boundary between 
survival and destruction for real stars is likely to be slightly 
different than what is predicted here, as the central densi- 
ties of stars on the MS depend on rotation, metallicity, and 
age (Maeder 1974; Wagner 1974). Notably, our own Sun 



has a central density approximately twice that of the standard 
7 = 4/3 polytrope used to model it. This may allow the cores 
of somewhat evolved MS stars to survive for slightly larger 
values of (3, although their gravitational influence is likely 
small as the helium-enriched cores of evolved MS stars are 
no larger than 10% of the star's mass (Schonberg & Chan- 
drasekhar 1942). 



3.2. Characteristic features ofMit) 

Figure 5 shows the family of M(t) curves as a function of 
f3 for both 7 = 4/3 and 7 = 5/3. Immediately evident is the 
strong dependence between M pea k and j3 for (3 < /?d, and the 
similarity of the M(t) curve family for (3 > (3^. The result 
that deeper encounters do not produce more rapid flares is 
in direct conflict with the analytical prescription presented in 
Lodato et al. (2009) (hereafter LKP), in which the binding en- 
ergy dM/dE is equivalent to the spread in mass over distance 
(modulo a constant), dM/dx, at pericenter. In this model 
(hereafter referred to as the "freezing model"), the binding 
energy is given by 

E = GM h x/r 2 p , (4) 

and thus deeper encounters always result in faster-peaking 
transients. Because the binding energy E oc r~ 2 , the scaling 

between (3 and f pea k is expected to be f pea k oc /3 3 (Ulmer 1999). 

We definitively find that this is not the case, as the 
two separate functional forms of the parametric pair 
[fpeak(/?),^peak(/?)] indicate a separate set of assumptions are 
appropriate for the two cases j3 < /3<j and j3 > f3i, neither of 
which match the functional form advocated by LKP (Fig- 
ure 5, triangles). For encounters in which j3 < /3d, fpeak and 
Wpeak are approximately related to one another by a power 
law, with the best fit model having M pea k oc f p( J a £ for 7 = 4/3 
and M pea k oc fp ea ° k 5 for 7 = 5/3. The steepness of this rela- 
tion means that the difference in f pea k is only a few tenths of 
a dex between an event in which 1O~ 4 M is lost and a full 
disruption. For j3 > /3 d , the trend between f pea k and M pea k re- 
verses for increasing j3, with deep encounters resulting in both 
slightly longer duration flares and slightly lower typical accre- 
tion rates. 

For fully-disruptive encounters, we find that M(t) varies lit- 
tle with increasing (3. An assumption of the freezing model 
is that the distance at which the dynamics of the debris can 
be described by Kepler's laws is when the star is at pericen- 
ter. In fact, the star's self-gravity becomes unimportant be- 
fore the star comes this close to the black hole for encounters 
where f3 > j3^. This suggests that the binding energy distri- 
bution of the material should be determined shortly after the 
star crosses the full disruption radius r^ = r t //3d, and not at its 
closest approach, unless the encounter is grazing enough such 
that r p < r d . 

This can be understood by considering the local reaction 
time of each layer of the star's structure as compared to the 
passage timescale. The dynamical timescale for a particular 
layer is r d y n — y/l/Gp x , which is approximately equal to the 
time between when the star is at a distance where the tidal 
force is capable of removing that layer and the time of peri- 
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FIG. 5. — Fallback accretion rate M(t) onto a 1O 6 M0 black hole from the disruption of a 1 Mq star as a function of 7 and 0. The colored curves in the left 
two panels show M(f), with the color of each curve corresponding to the color coding scheme presented in Figures 1 and 2. The dashed portions of each curve 
are extrapolations based on the slope of the M(t) immediately prior to the extrapolated region, which accounts for the fact that the late-time slope can only be 
determined exactly if the simulations are run for a prohibitive amount of time (see Figure 10). The dotted line shows the Eddington limit for a 10 6 Mq black 
hole assuming an accretion efficiency e = 0. 1 . The open triangles connected by the gray dashed line show the peak fallback rate M pc;u \; and time of peak t™.^ as 
predicted by the energy-freezing model, in which the period of the return of materials scales as 1 (Evans & Kochanek 1989; Ulmer 1999; Lodato et al. 2009). 
The open circles connected by the black line show the fits to and ? pea i; as given by equations Al and A2 respectively. The right two panels shows the same 

data as the left two panels, with the filled regions showing the range of M (t) curves resulting from full disruptions (red) and from disruptions in which the star 
survives (gray). 

forces continuously evolves over the encounter, the actual ra- 
dius at which mass is removed can be larger or smaller than 
r e g, and thus the relationship between E and x is more com- 
plicated than outlined here. 

Additionally, while the binding energy is effectively frozen- 
in once the star crosses r<j, the assumption that the orbital en- 
ergy can be reliably recorded at this point is only valid if the 
pressure gradient that develops within the star during maxi- 
mum compression is not large enough to affect dM/dE. As 
shown in Carter & Luminet (1983), the pressure component 
of the Lagrangian does build significantly shortly after peri- 
center, and eventually dominates the tidal component for suf- 
ficiently deep encounters. However, while this build-up can 
lead to the production of shocks whose breakouts may be ob- 
servable as short X-ray transients (Guillochon et al. 2009), we 
find that the gradient of pressure within the orbital plane pri- 
marily acts to redistribute the most highly-bound material for 
(t < ?peakX an d not the material that determines the behavior 
of the decay phase (see Figure 3 of Guillochon et al. and Fig- 
ure 5 of LKP). The tangible effect of this pressure build-up 
on the shape of M(t) is the spreading of some material that 
would have otherwise accreted at f pea k to more highly-bound 
orbits, thus reducing the rate of accretion at peak, shifting f pea k 
to later times, and leading to an increased feeding rate at early 
times. 

This behavior is analogous with what is found for binary 
star disruptions, in which the change in orbital energy AE OI b 
of the stars is independent of the impact parameter for suf- 
ficiently deep encounters (Sari et al. 2010). In Figure 6, we 
compare AE^ calculated by Sari et al. for binary disruptions 
to the mass-averaged spread in the binding energy (E) of the 
material that becomes bound to the black hole. As is found in 
binary disruption calculations, the change in energy initially 
increases with increasing (3, then a transition point is reached 
where the binary's gravity (or star's gravity, in our case) no 
longer affects the dynamics, and finally the change in energy 
approaches a constant. A single star disruption and the dis- 
ruption of a binary system are conceptually quite similar. A 
full disruption is analogous to an equal-mass binary disrup- 
tion with separation distance a ~ R*, where the mass of each 
"star" is equal to the mass liberated from each Lagrange point, 



FIG. 6. — Average spread of matter post-disruption as compared to the 
change in orbital energy for binary disruptions (Sari et al. 2010). The solid 
curves show (E) (the mass-averaged binding energy of the bound debris 
post-disruption) for both 7 = 4/3 (light blue, open circles) and 7 = 5/3 (or- 
ange, filled circles) stars, whereas the dashed curves show the change in 
orbital energy A£ or i, for prograde binary encounters, where we have pre- 
sumed that each star has mass Mq . The impact parameters f}$ beyond which 
stars are considered to be destroyed being denoted by the colored dot-dashed 
lines. The red dashed curve shows the maximum change in orbital energy 
A£ ol -b.max at a particular f3, whereas the blue dashed curve shows E OI \, av- 
eraged over binary phase. The binary disruption energies are scaled by 
(GM t /a)(M h /M Q y/ 3 , where a is the initial binary separation, whereas the 
stellar disruption curves are scaled by (GM* /R„)(Mh/M0)' / ' 3 . 



center, 



Ttidal = r t ,x/v t ,x — 



GM h V GM * 



(5) 
(6) 



where the subscript x refers to quantities defined by the mass 
interior to x. Thus, regardless of the distance at which the tidal 
force begins to dominate the self-gravitational force, material 
is removed from the star at or near the full disruption radius rj. 
This means that the effective radius that should be used in the 
denominator of equation (4) is r e g = max(rd,r p ). However, as 
the degree of balance between the tidal and self-gravitational 
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FIG. 7. — The left panel shows the power-law index n, with orange corresponding to 7 = 4/3 stars and light blue corresponding to 7 = 5/3 stars. The filled 
regions show the convex hull of n over all values of 0\ the lightly-weighted curves within the filled regions represent individual simulations for specific values 
of 0. The horizontal dotted and dashed black curves show n = (i.e. the peak of the accretion rate) and n = —5/3 (the canonical value for constant dM/dE), 
respectively. The dashed colored curves are produced using the analytical formulae of Lodato et al. (2009) for = 1 encounters for both values of 7. For 1 Mq 
stars, the analytical formulae predict a faster rise to peak for 7 = 4/3 (orange-bordered point) than for 7 = 5/3 (blue-bordered point), and fail to reproduce the 
steeper power law index that is found shortly after peak for 7 = 4/3. The open colored circles show the numerical results of Lodato et al. (2009) for = 1, where 
we include their 7=1.4 case (dark green) in addition to 7 = 5/3 and note that their simulations set V = 7. The right panel shows the asymptotic power-law index 
fioo as a function of 0, with the color coding scheme identical to that of Figures 1 and 2, where the filled circles show «oo for 7 = 4/3, and the open circles for 
7 = 5/3. A best fit for both values of 7 are shown by the solid colored curves, and the impact parameters 0$ beyond which stars are considered to be destroyed 
are denoted by the colored dot-dashed lines. 



M\ = Mi = M*/2. A partial disruption is analogous to an un- 
equal mass "trinary" system, in which the three masses corre- 
spond to the surviving self-bound core with mass - AM, 
and the bound/unbound debris streams with mass AM/2, all 
with initial separation a ~ R*. As the results presented in Sari 
et al. are independent of mass, the normalized AE for dis- 
ruptions still map closely to those seen for binary disruptions, 
despite the variance in mass of the two (or three) interacting 
objects. 

A caveat in our comparison to binary systems is that bina- 
ries can have an arbitrary orientation upon arrival at percen- 
ter. This leads to an increase in the potential maximum en- 
ergy change, and the f3 value at which it occurs. This comes 
as the result of stars in a binary being able to come arbitrarily 
close to one another during an encounter for favorable binary 
phases at pericenter (A£'orb,max m Figure 6), which permits 
them to interact gravitationally for longer. In effect, the bi- 
nary can become "denser," decreasing the size of its effective 
tidal radius. For a stellar disruption in which the star is not 
initially rotating, the mass interior to a given radius cannot 
increase in the same way, and thus its self-gravity ceases to 
be important interior to its original tidal radius. This results 
in the near-constant spread in energy as described above, and 
is visually evident from simulation snapshots (Figures 1 and 
2), in which the debris distributions are almost identical for 
^ /3d- We speculate that for rapidly rotating stars that this 
same effect that can yield large AE for certain binary phases 
may also apply, as stellar rotation can permit stars to penetrate 
more deeply before dM/dE is set (Stone et al. 2012). 

Figure 7 shows both the power-law slope 71(f) over the full 
M(t) curve (left panel), and the asymptotic power-law slopes 
"00 (right panel), as produced by our disruption simulations. 
The behavior of the curves is more complicated than what is 
implied by the freezing model, in which n > -5/3 for all t. 
The qualitative behavior of the M(t) curves can be character- 
ized by three phases: A rise phase, in which n > 0, a drop 
phase, in which n < (and potentially even < -5/3), and an 
asymptotic phase, in which = n(t — > 00) assumes a con- 
stant value. The rise phase is somewhat similar to what is 



predicted by the freezing model, although the evolution of n is 
somewhat more rapid within the simulations. The drop phase 
exhibits particularly steep downward slopes for 7 = 4/3 stars, 
with shortly after peak, but n for 7 = 5/3 stars is closer 

to the predicted asymptotic value. Despite the disagreement 
between our simulations and the analytical model presented in 
LKP, we do find we reproduce the simulation results of LKP 
for 7 = 5/3 stars at the same /3 (Figures 7 and 13). 

The discrepancy between the simulation results and the pre- 
diction of the freezing model can be understood by consid- 
ering what material becomes bound to the black hole in the 
case that the star is not completely destroyed. The bind- 
ing energy E of material a distance x from the center of 
the star is cx x/ max(rd,r p ) 2 oc xf3 2 (assuming x -C r p ). The 
value of x that corresponds to the material that determines the 
asymptotic behavior of M{t) can be estimated by consider- 
ing the deepest point within the star during the encounter in 
which tidal forces are capable of overcoming the star's self- 
gravitational force, the exact functional form of which is de- 
pendent upon the hydrodynamical response of the star during 
the encounter. While this functional form can only truly be de- 
termined through hydrodynamical simulation, it is clear that 
the effective x must decrease with increasing 0, as the tidal 
forces remove an ever-increasing fraction of the star's mass. 
This implies that the scaling between E and ft must be weaker 
than f3 2 , and thus f pea k should show less evolution for progres- 
sively deeper, but not completely disruptive encounters. 

The asymptotic phase exhibits a more complicated behavior 
that depends on (3, and shows four distinct behaviors depend- 
ing on the depth of the encounter for both 7 = 4/3 and 7 = 5/3 
stars. For extremely grazing encounters in which a small frac- 
tion of the star's mass is lost, ~ -5/3. In these encounters, 
all of the mass is removed near pericenter, resulting in an en- 
ergy spread that only depends on x, the distance to the star's 
center of mass (in agreement with the freezing model). 

When a significant fraction of the star's mass is removed 
in an encounter, steepens to values as large as ~ -2.2. 
This behavior arises from the influence of the star's core (Fig- 
ure 8). As the outermost layers of the star are removed prior 
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FIG. 8. — Cartoon showing the gravitational effect of a surviving core on 
the dynamics of material that is removed from the star during a partially- 
disruptive encounter. The inset diagram in the upper left demonstrates how 
the restoring force provided by a surviving core can alter the structure of the 
outer layers. For encounters in which the core plays little role, the binding 
energy to the black hole E scales linearly with x (black solid curve), the dis- 
tance from the star's center of mass (Lodato et al. 2009). If a core survives the 
encounter, its gravity prevents material from moving as quickly away from 
the star, resulting in a weaker relationship between E and x (orange dashed 
curve). This consequently results in a more-steeply declining M(t). If the 
core itself is close to destruction, its gravitational influence is minimal, and 
the canonical M(t) <x f~ 5 ' 3 decay law is recovered. 

to pericenter, the core is able to partially counter the black 
hole's tidal force, keeping material closer to the star's core, 
and thus reducing the effective x at which the material is no 
longer strongly affected by the core's gravity. This results in 
a sub-linear relationship between E and x. As E oc f" 2 / 3 , and 
E oc x' n , where m < 1, the resulting asymptotic power-law is 



where = -5 /3 is recovered for the standard linear relation- 
ship. If m < 0, this implies that E actually decreases with x, 
and thus the most bound material would initially lie interior 
to the least bound. As the most bound material would have to 
cross beyond the least bound, we expect that any m < rela- 
tionship would be quickly flattened to at least m = by pres- 
sure gradients, resulting in a limit on the asymptotic slope of 
«oo > -7/3. 

For disruptions that are just deep enough to destroy the star 
(i.e. /3 = Pa), n^ can be somewhat less steep than -5/3. This 
implies that m > 1 ; the relationship between E and x is super- 
linear. For this borderline case only, the release of material 
that eventually composes the decay tail of M(t) is moderated 
by the slow shrinkage of the stellar core, which is not fully 
destroyed until after the star has passed pericenter. For these 
encounters, r and x are somewhat dependent, with the material 
being released at the smallest x being launched at large r, and 
thus the quantity x/r 2 can be oc x >l . 

Finally, for deep encounters, again seems to be consis- 
tent with -5/3. Unlike the borderline case where a core per- 
sists long after pericenter, here the core is rapidly destroyed, 
and the energy is again set at a fixed r. As there is no core 
to resist the tidal force, the energy spread is simply given by 
the spread in potential energy across the star, a la the freezing 
model. 



3.3. The Influence of Stellar Structure 

A seemingly counter-intuitive result in the context of the 
freezing model is the fact that less-centrally concentrated 
stars, which have more mass at larger radii, result in transient 
events that peak at later times than their more centrally con- 
centrated brethren, even for events that yield the same AM. 
As shown in Section 3.1, the distance at which total disrup- 
tion occurs is significantly deeper for centrally concentrated 
stars, which may explain some of the discrepancy. Consider 
what happens to a star in the approach to pericenter for two 
extreme cases: A case in which most of the star's mass is 
concentrated at its center, and a case in which the star has 
near-constant density. In the centrally-concentrated example, 
the star's outer layers will find that their dynamics are partly 
determined by the tidal force at early times, but also partly 
determined by the core, which remains initially undisturbed 
(Figure 8). The influence of any surviving core on the dy- 
namics of the matter can thus affect the final binding energy 
E. 

The left panel of Figure 9 shows that one of the fundamen- 
tal assumptions of the freezing model, that the binding energy 
E and the distance from the star's core x are linearly related, 
is not correct, and matter that contributes to a particular E 
is drawn from a range in x that spans nearly the entire star. 
In Figure 9 we compare two disruptions with nearly identical 
AM for 7 = 4/3 and 7 = 5/3. The left panel shows the time 
evolution of the material that determines the peak of M(t), 
with the contour colors corresponding to the same times af- 
ter pericenter. The filled contours show the regions that con- 
tribute to Mpeak within each snapshot. The expectation un- 
der the freezing approximation would be that all mass that 
possesses a given energy E comes from a cylindrical cross- 
section of the star, but as illustrated in Figure 9 the geometry 
of the debris that contributes to M pea k is clearly not cylindrical. 
The right panel shows the M(t) that correspond to these snap- 
shots. For the earliest snapshots (light blue-filled contours), 
the material in the 7 = 5/3 simulation appears to have a head- 
start over the centrally-concentrated case, despite the fact that 
the 7 = 4/3 encounter is 50% deeper ((3 = 1.0 for 7 = 5/3 
versus (3 = 0.65 for 7 = 4/3). However, as the encounter pro- 
gresses, the peak of M(t) for 7 = 4/3 moves to progressively 
earlier times relative to 7 = 5/3, eventually settling to a value 
that results in a faster transient. 

This implies that the binding energy of these layers relative 
to the black hole should not be recorded assuming the star 
has preserved its original spherical shape and size. While less 
material is positioned near the black hole when comparing 
the centrally-concentrated case to the constant density case, 
the core of the star continues to interact with the debris, and 
effectively "carries" material to larger E before E has been 
fixed. For constant density stars, the tidal force and the self- 
gravitational force scale to the same power in x, and the core 
is disrupted at approximately the same time as the outer lay- 
ers, which is consequently why these stars are destroyed at a 
distance that more closely matches the classical Roche result. 
In this case, the assumption that E can be determined by con- 
sidering the star's original size and shape is more appropriate, 
as little stellar material is carried closer to the black hole, as 
is found in the centrally-concentrated case. 



3.4. Long-term evolution ofMit) 



10 



Guillochon, Ramirez-Ruiz 





-0.8 -0.6 

1 Log 10 f (yr) 

FIG. 9. — The left panel shows a collage from two simulations of the regions that contain the mass that contributes to M ptil ^(t), super-imposed over the 
density distribution of the each snapshot. The fill and outline color denotes the time of the snapshot after pericenter, with blue corresponding to t = f p , and red 
corresponding to t = t p + 2 X 10 4 . Snapshots are shown from two different disruption simulations that have similar values of AM, with light blue showing the 
disruption of a 7 = 5/3 star for (5 = 0.65, and orange showing the disruption of a 7 = 4/3 star for j3 = 1.0. Each fill region shows the material that contributes to 
the part of M(t) that is within 90% of M pe ak(0- The right panel shows the values of M(t) derived from these two simulations (one curve per snapshot), with the 
colors of each curve corresponding to the color of the fill regions of each snapshot. The solid lines correspond to 7 = 5/3, and the dashed lines corresponding to 
7 = 4/3, with the thick light blue and orange lines being fitted to M pca ] t (f). 



While the peak of the accretion rate is determined within 
tens of stellar dynamical timescales, the tail of M{t) contin- 
ues to evolve for hundreds of dynamical timescales. Most 
notably, a large "cavity" in M(t) is present for accretion times 
t > (Figure 10), which gradually fills from left to right as 
additional material in the tidal debris tails satisfies the simple 
energy criteria applied to determine if matter is bound to the 
black hole. An example of the long-term evolution of M(t) 
arising from this interaction is shown in Figure 10, where 
M(t) is not determined for r > 10 years until 550 dynami- 
cal timescales after the disruption. As the star recedes from 
the black hole, the evolution of M(t) slows (as indicated by 
the decreasing space between the curves in Figure 10), imply- 
ing that progressively longer simulation times are required to 
determine M(t) much beyond f pea k. 

The cavity arises from the exclusion of material within the 
debris stream that remains bound to the stellar core after the 
encounter, with the evolution coming as a result of the contin- 
ued interaction between either the stellar core (if the star sur- 
vived the encounter) or a mildly self-gravitating debris stream 
(if it did not) and the black hole. An examination of the pres- 
sure of the debris tails reveals that the debris is free-streaming, 
even in the vicinity of the Hill sphere. In other words, the 
pressure gradients present within the stream are small enough 
to be incapable of modifying the material's trajectory. Thus, 
the interaction between the stream, black hole, and surviving 



core is purely gravitational. 

In Figure 1 1 we show dM/dE for disruptions of a 7 = 4/3 
star for three values of f3. As material that is considered to 
be bound to the star (yellow) crosses the time-dependent Hill 
sphere, it becomes bound to the black hole (green). We find 
that there is always some mass in the vicinity of the time- 
dependent Hill sphere «h(0 (cyan) for encounters in which 
a core survives, whereas full disruptions do not show this 
behavior, mostly because the self-bound mass shrinks dras- 
tically as progressively less material satisfies the criteria for 
being self-bound. 

As the star moves away from the black hole, the Hill radius 
grows, but by a rate that is a factor (M^/Mh) 1 / 3 smaller than 
the rate at which the star recedes from the black hole. This 
implies that any material that retains a positive velocity rela- 
tive to the surviving core after the encounter has the potential 
to be removed from the star, even if it is technically bound 
to the core (i.e. v 2 < 2GM core /V) at an earlier epoch. Much 
of the observed evolution of dM/dE may be due to our def- 
inition for what is considered to be "bound" to the surviving 
core after the encounter (Section 2). Our calculation presumes 
that the energy budget of material with respect to the star is 
sufficient to determine what inevitably remains bound to the 
star; in reality the question of whether a given particle remains 
bound or not amounts to solving the restricted elliptical three- 
body problem, for which no closed-form solution exists. The 
Jacobi constant, which has a fixed value in the restricted cir- 
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FIG. 10. — Evolution of dM/dE (left panel) and M(t) (right panel) for the disruption of a star with polytropic index 7 = 5/3 and impact parameter f} = 0.55. 
The measured mass distribution in both plots is shown as a function of time relative to the time of pericenter, with the blue curves indicating early times and the 
red curves showing late times. The right-hand plot is overlaid on a series of dotted gray lines showing the fiducial M(t) <x r~ 5//3 evolution. While the early-time 
M(t) can be determined shortly after the disruption, the rate of fallback still evolves for t > 10 yr even after the simulation has been allowed to run for many 
hundreds of dynamical timescales. 
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FIG. 1 1 . — Distribution of mass as a function of binding energy E for three different simulations of a 7 = 4/3 disruption for j3 = 0.6, 1.4, and 2.5. Shown in 
each panel is a sequence of mass distributions in time, with the evolution in time progressing from top to bottom, where the yellow regions show all material that 
remains bound to the black hole, and the green regions show material bound to the black hole but not bound to the star. The cyan region shows material that is 
bound to the star, but whose semi-major axis larger than the distance defined by the surviving core's time-dependent Hill sphere, au(t). The green regions are 



separated by the gray dotted curves into sub-regions where the resulting accretion rate M(t) exceeds, from lowest to highest, 10 10 10 -4 , 10 3 , 10 1 , 10 
and 1 Mq/yt. The red arrows show the location within the material bound to the black hole which determines the peak accretion rate M^.^. 
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cular three-body problem and can be used to determine the 
zones within which a particle of a given initial position and 
velocity can occupy (Murray & Dermott 1999), is not con- 
stant once the orbit is non-circular (Hamilton & Burns 1992). 

However, while the energy balance approach may not be 
capable of immediately determining the mass that will even- 
tually become bound to the black hole, the distribution only 
remains uncertain for material that is accreted far beyond t p&: ±. 
In the limit that t — > 00, the distance of the star to the black 
hole increases as J 2 / 3 , and thus E for material leaving the Hill 
sphere is GM^a^/r 2 oc f" 2 / 3 . This explains the observed slow- 
ing of the evolution of M{t). 



The material that fills in the cavity assumes a distribution in 
E that is not entirely flat, resulting in a fallback rate that scales 
to a power slightly steeper than the canonical f" 5 / 3 . This en- 
ergy distribution is likely set near pericenter, where the pres- 
sure component is comparable to the tidal component. As 
the star recedes from the black hole, the pressure component 
of the force decreases more quickly than the tidal component 
(Kochanek 1994), and thus the debris is expected to evolve 
purely gravitationally. However, the conditions under which 
material is launched across the time-dependent Hill sphere 
may depend somewhat on the pressure gradient, no matter 
how small, as the net gravitational force is zero (Lubow & 
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Fig. 12.— Fits to M, 



peak 



and fp Ca i;, with the fits to the 7 = 4/3 models being shown by the solid colored circles, and fits to the 7 = 5/3 models being shown by 
the open colored circles. The color coding matches that of Figures 1 and 2, with the impact parameters 04 beyond which stars are considered to be destroyed 
being denoted by the colored dot-dashed lines. 



Shu 1975). 



4. DISCUSSION 



The results of our extensive parameter study produce a 
number of unexpected trends as compared to the predictions 
presented by previous work. In the previous section we at- 
tempted to explain the observed scalings, and how these fea- 
tures arise as a result of the interaction between the black hole 
and a potentially surviving stellar core. In what follows, we 
explain how these newly discovered features can be used to 
constrain the type of disrupted star and how it was disrupted. 

4.1. Can 7 and f3 be determined a posteriori? 

Given our predicted M(t), is it possible to determine either 
the stellar structure or the impact parameter from the light 
curve produced by a tidal disruption event? The conversion 
efficiency between mass accreted by the black hole and the 
light emitted is somewhat uncertain, and depends on factors 
such as the black hole's spin and how the accretion rate com- 
pares to M Edd (Ulmer 1999; Beloborodov 1999; Strubbe & 
Quataert 2009; Lodato & Rossi 2010). However, as the effi- 
ciency cannot be larger than unity, and as flares are typically 
observed in the decay phase, we can only place lower limits 
on the amount of mass accreted by a black hole to produce a 
given flare (Gezari et al. 2008). Thus, at the very least, our 
predicted AM (Figure 4) can be used to exclude events for j3 
less than some critical value, given the mass of the star. 

Figures 7 and 12 present four additional quantities that en- 
able us to classify tidal disruptions based on the properties of 
observed tidal disruption flares. Two of these quantities, M pea k 
and fp ea k, are only available to us for flares in which the peak 
of the accretion rate is clearly observed (Gezari et al. 2012), 
but both n(t) and n^ are measurable for flares that are ob- 
served long after peak (Komossa & Greiner 1999; Komossa 
et al. 2004; Gezari et al. 2006, 2008; Cappelluti et al. 2009; 
van Velzen et al. 201 1; Cenko et al. 2012). If the mass of the 
black hole is known with some certainty, one may be able to 
infer both and j3 by simply measuring and f pea k and 
comparing to our resultant M(f), which at first glance appear 
to form distinct sequences for 7 = 4/3 and 7 = 5/3 stars (Fig- 
ure 13, left panel). However, this is only true assuming that 
centrally concentrated stars have the same mass and radius as 
stars of near constant density. The transition from stars that 
are well-modeled by a 7 = 4/3 poly trope to a 7 = 5/3 poly- 
trope is also accompanied by a decrease in radius such that all 



stars with mass O.25M < < M Q have the same central 
density (Kippenhahn & Weigert 1990). Adjusting the radii 
and mass of our 7 = 5/3 models to the mass and radius of 
a 0.25M© star (Tout et al. 1996), we find that the sequence 
of M(t) functions for 1.0M Q and O.25M stars lie on top of 
each other (Figure 13, right panel), making the determination 
of the disrupted mass of a star somewhat degenerate with its 
structure. 

This motivates us to look for other features of M(t) that 
may uniquely identify either 7 or /?. If we consider the power- 
law of the rate of decline n after peak, we find that there is a 
distinguishing feature between 7 = 5/3 and 7 = 4/3 models at 
~ 0.5 dex after fp ea ^. Whereas 7 = 5/3 stars quickly converge 
to n ~ -5/3, 7 = 4/3 models show a characteristic drop, with n 
being as large as -4 for some encounters (Figure 7, left panel). 
This feature is most prominent for intermediate /3 in which 
~ 50% of the star's mass is removed during the encounter, and 
represents the strong influence of the dense stellar core, which 
acts to drag material deeper within the black hole's potential 
before tidal forces are capable of removing it. 

In addition to being more centrally-concentrated to begin 
with, an additional component that likely contributes to this 
observed drop is the adiabatic response of the surviving core. 
For 7 = 5/3 stars, the removal of mass results in the infla- 
tion of the star, whereas 7 = 4/3 exhibit the opposite behav- 
ior, shrinking dramatically in response to the loss of mass 
(Hjellming & Webbink 1987). This enhances the core's in- 
fluence during the encounter in the phase where the core's 
mass is changing, slowing the reduction in the core's effective 
gravity, and thus pulling even more matter to higher binding 
energies. The recently observed flare PSl-lOjh presented in 
Gezari et al. (2012) shows a clear drop in the accretion rate 
with respect to the canonical f" 5 / 3 decline rate expected from 
the freezing model. In the freezing model, it is impossible to 
produce a decline feature steeper than t~ 5 / 3 within any part 
of M(t), as we explained in Section 3.2. As many tidal dis- 
ruption flares may show this characteristic drop in M(t), a 
clearly-resolved peak can be used to compare to the subse- 
quent decay phase for a precise determination of n(t). 

For events in which the peak is not clearly observed, and 
for which the signal-to-noise is too small to permit an accu- 
rate determination of n(f), the asymptotic slope n^ of M(t) 
can still provide additional information about the star that was 
disrupted. As shown in the right panel of Figure 7, n^ can 
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FIG. 13. — Comparison of the families of M(t) curves for 7 = 4/3 (orange) and 7 = 5/3 (light blue). The left panel shows M(t) derived from the simulations 
presented here as solid lines, assuming that the stars for both 7 have mass M = IMq. The open circles show the numerical results of (Lodato et al. 2009) for 
= 1 and 7 = 1.4 (dark green) and 7 = 5/3 (light blue). The right panel shows the shift in M(t) (black arrows) if the mass and radius of star that is expected to 
have a structure described by 7 = 5/3 is taken into account (Tout et al. 1996). 



be used to distinguish between partial and full disruptions. 
The fact that assumes values that are significantly steeper 
than -5/3 may indicate that additional tidal disruption flares 
have been found observationally, but subsequently discarded 
and/or ignored due to the mismatch between the measured n 
and —5/3 (van Velzen et al. 2011). This implies that some 
supernovae that have been observed at the centers of galaxies 
may in fact be misidentified partial tidal disruptions. 

4.2. Future work 

As found in previous work (Faber et al. 2005; Guillochon 
et al. 2011), there is a change in surviving star's orbital en- 
ergy after the encounter, with the change in energy being 
comparable to the star's initial self-binding energy. This 
change in energy, combined with the star's initial orbital en- 
ergy, leads to a shift in the entire dM jdE distribution, which 
can affect the fallback of material for E ~ A£ or b, or for t > 

M h Rl /2 G-^ 2 M~ 3/2 ~ 100 years, given that A£ orb ~ GM»//?». 
As the star's initial orbital energy may not be zero and can be 
comparable to AE itself, and thus the final binding energy 
of the star depends on its initial orbital energy, we have pre- 
sented our dM/dE and M(t) curves with this change in en- 
ergy removed. As a result, our plots show the fallback rate 
that would be expected if the final star were to remain on a 
parabolic trajectory, as our initial conditions assume. While 
these kicks that are typically of the order of star's own escape 
velocity may be important in determining the further fate of 
the star and whether it will suffer additional disruptions, they 
are not expected to affect the first century of a flare's evolu- 
tion, of which only the first few years are accessible to cur- 
rently available transient surveys. 

Even if M(t) is directly related to the properties of the star 
being disrupted, the luminosity of the accretion disk L may 
not directly follow M{t). The primary factors that affect the 
link between M{t) and the bolometric L are the viscous evo- 
lution of the disk and the size of the disk (Ramirez-Ruiz & 
Rosswog 2009), although other processes may strongly af- 
fect the amount of light observed in a single band, especially 
in the optical/UV where dust extinction can play a vital role. 
Disk viscosity can only affect L for t < r v i sc , in which its pri- 
mary affect is to delay emission at early times. However, once 
t > T v i sc , L is expected to track M(t) closely. As the material 



is delivered to the disk at r ~ r p , the ratio of r v j sc to f pea k is 
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where Z? 7 is a fitted parameter derived from our simulations, 
(B 7 ^0.1 for most /3, see Appendix A) and a is the param- 
eterized a-disk scaling coefficient, where we have taken the 
scale-height ratio h/r=\. If T v i sc /f pea k > 1, the accretion is 
spread over longer timescales, resulting in a power-law decay 
index n = -1.1 (Cannizzo et al. 1990). This may affect the 
light curve shape in the earliest phases of the fallback (prior 
to peak) where t <C r v i sc , and thus the early evolution of L(t) 
may not follow the functional forms of M(t) presented here 
for t <C fp ea k. But as observations of tidal disruption flares 
in the decay phase seem to be consistent with the canonical 
n = -5/3 decay law, it is clear that M(t) and L(t) must be 
closely coupled on year-long timescales. 

An ingredient that the set of simulations presented in this 
paper do not include is the inclusion of general relativistic 
effects, which become important for very deeply penetrating 
encounters. Qualitatively, for both spinning and non-spinning 
black holes, general relativity is expected to result in more 
mass loss and a spreading of mass in dM/dE as compared 
to Newtonian encounters, as its primary effect is to bend the 
star's path such that it spends a larger fraction of time near the 
black hole where tidal forces are important (Luminet & Marck 
1985; Kobayashi et al. 2004). The only numerical prove- 
nance for how the metric may affect the feeding rate comes 
from low-resolution simulations performed by Laguna et al. 
(1993), which find a slight increase in M pea k for increasing [3, 
but much less than the predicted /3 3 scaling. If the black hole 
has non-zero spin, the resulting dM/dE depends on the ori- 
entation of the star's angular momentum vector as compared 
to the black hole's spin vector (Haas et al. 2012). 

A spinning black hole permits deeper encounters that 
don't result in the star being immediately swallowed (Kesden 
2012b), provided that the two angular momentum vectors are 
aligned, and also should affect the final binding energy dis- 
tribution, with co- and counter-rotational encounters resulting 



14 



Guillochon, Ramirez-Ruiz 



in smaller and larger AM, respectively (Diener et al. 1997; 
Ivanov & Chernyakova 2006; Kesden 2012a). However, as 
the fraction of disruptions in which non-Newtonian metrics 
can affect the dynamics is ~ r s /r t , which is ~ 5% for a 10 6 M Q 
black hole and ~ 20% for a 10 7 M Q black hole, the majority 
of tidal disruption events are well-represented by a Newtonian 
approximation to the black hole's gravity. 

Lastly, the absence of hydrogen in spectra taken of the tidal 
disruption event PSl-lOjh (Gezari et al. 2012) strongly sug- 
gests that the disruption of stars that are not on the MS may 
contribute significantly to the overall rate of tidal disruption. 
As we show, the structure of the star that is disrupted is clearly 
imprinted upon M(f), providing valuable additional informa- 
tion that can be used to distinguish between candidate disrup- 
tion victims. We explore the disruption of post-MS stars in a 
companion paper using a method similar to what is presented 
here (MacLeod et al. 2012). 

The discovery of flaring black hole candidates in nearby 
galaxies will continue to elucidate the demography of the 
AGN population (De Colle et al. 2012). Whereas AGN are 
supplied by a steady stream of fuel for hundreds or even thou- 
sands of years, tidal disruptions offer a unique opportunity to 
study a single black hole under a set of conditions that change 
over a range of timescales. There are, of course, rapidly vary- 
ing stellar-mass black hole candidates in X-ray binaries within 
our own Galaxy. But for SMBHs, tidal disruption events offer 



the firmest hope of studying the evolution of their accretion 
disks for a wide range of mass accretion rates and feeding 
timescales. The simulations and resultant M(f) curves pre- 
sented here are crucial for determining the properties of the 
black hole itself, as an incomplete model of a stellar disrup- 
tion can result in much uncertainty in how the black hole con- 
verts matter into light. For a disruption with a well-resolved 
light curve, our models permit a significant reduction of the 
number of potential combinations of star and black hole prop- 
erties, enabling a better characterization of SMBHs and the 
dense stellar clusters that surround them. 

We have benefited from many useful discussions with S. 
Gezari, M. Macleod, C. Miller, S. Liu, R. O'Leary, F. Ra- 
sio, M. Rees and C. Thompson. We thank the anony- 
mous referee for constructive corrections and suggestions. 
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the DOE-supported ASCI/ Alliance Center for Astrophysical 
Thermonuclear Flashes at the University of Chicago. Com- 
putations were performed on the UCSC Pleiades and Laozi 
computer clusters, and the NASA Pleiades computer cluster. 
We acknowledge support from the David and Lucille Packard 
Foundation, NSF grants PHY-0503584 and ST-0847563 and 
the NASA Earth and Space Science Fellowship (JG). 



APPENDIX 



FITTING PARAMETERS 

For convenience, we have calculated fitting parameters for four characteristic quantities: The peak accretion rate M pea k, the 
time of peak accretion f pea k, the amount of mass lost by the star AM, and the asymptotic decay power-law index n^. These 
parameters can be used to constrain observed tidal disruption events based on measurable characteristics of their light curves (see 
Section 4.1): 
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In these expressions are four functions of j3 alone: A 7 , B~ l , C 7 , and D~ ( . The forms of these functions are derived by fitting ratio- 
nal functions to the outputs produced by the numerical simulations presented in this paper. These functions are derived separately 
for two polytropic 7, 7 = 4/3 and 7 = 5/3, which are appropriate for high- and low-mass main sequence stars, respectively. 
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